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Abstract. We present the first "combined stellar struc- 
it«je and atmosphere models" (CoStar) for massive stars, 
^^lich consistently treat the entire mass loosing star from 
©te center out to the asymptotic wind velocity. The mod- 
T et^ use up-to-date input physics and state-of-the-art tech- 
toQues to model both the stellar interior and the spheri- 
£3lly expanding non-LTE atmosphere including line blan- 
•^^ting. Our models thus yield consistent predictions re- 
^^rding not only the basic stellar parameters, including 
^s^undances, but also theoretical spectra along evolution- 
ary tracks. On the same ground they allow us to study the 
^jfluence of stellar winds on evolutionary models. 
CO In this first paper, we present our method and investi- 
T gcl!te the wind properties and the interior evolution on the 
^H^ain sequence (MS) at solar metallicity. 

The wind momentum and energy deposition associ- 
^t>ed with the MS evolution is given and the adopted wind 
"pr operties are discussed. From our atmosphere calcula- 
tion s, which include the effect of multiple scattering and 
lfae overlap, we also derive theoretical estimates of mass 
fgss driven by radiation pressure. These values are corn- 
ered with the predictions from recent wind models of 
ttfis Munich group (Pauldrach et al. 1990, 1994, Puis et 
al. 1995). We find an overall agreement with most of their 
results. In addition, our models are better in reproducing 
the strong wind momentum rates observed in supergiants 
than those of Puis et al. 

A comparison between boundary conditions given by 
the conventional plane parallel and the new spherically ex- 
panding atmosphere approach is made. For the MS evolu- 
tion the evolutionary tracks and the interior evolution are 
found to be basically unchanged by the new treatment 
of the outer layers. However, for stars close to the Ed- 
dington limit, a small uncertainty in the behaviour of the 
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deep atmosphere is found which might marginally affect 
the evolution. Given the small spherical extension of the 
continuum forming layers in the considered evolutionary 
phases, the predicted stellar parameters differ negligibly 
from those obtained using plane parallel atmospheres. 

Key words: Stars: atmospheres - early-type - evolution 
- fundamental parameters - Hertzsprung-Russel (HR) di- 
agram - mass-loss 



1. Introduction 

This is the first paper in a series dealing with "com- 
bined stellar structure and atmosphere models" . In the 
present publication we introduce our method and study 
wind properties and the interior evolution on the main se- 
quence. The second paper in this series will deal with the 
spectral evolution of the models presented here (Schaerer 
et al. 1995, hereafter Paper II). A first study covering the 
Wolf-Rayet phases is presented in Schaerer (1995a, b). The 
CoStar models use up-to-date input physics and are based 
on state-of-the-art techniques for both the interior and at- 
mosphere modeling. The latter, in particular, includes the 
effects of line blanketing in the non-LTE expanding at- 
mosphere, as presented by Schaerer k Schmutz (1994a, b, 
henceforth SS94ab). 

Mass loss through stellar winds is (i) the determin- 
ing process for the evolution of the most massive stars 
(Mi ^ 2OM ; cf. Chiosi k Maeder 1986, Maeder k Conti 
1994), and is also responsible for (n) profoundly shaping 
the emergent spectral energy distribution of these stars 
(e.g. Kudritzki k Hummer 1990, Schmutz et al. 1992). 
Since the mass outflow is important in the outermost lay- 
ers, where radiation decouples from matter, stellar winds 



cal depth r < can be derived, in particular also the 
photospheric radius R T=2 / 3 at tr oss = 2/3. The boundary 
conditions are given by the velocity 



v(M r = M) = v* 



(2) 



which, given the atmospheric structure (see below), is de- 
termined by the above requirement for the optical depth. 
For the interior below expansion is neglected (but see 
Schaerer 1995a, b for Wolf-Rayet phases). From the conti- 
nuity equation we obtain the boundary condition for the 
density i.e. 



p(M r = M) = 



M 



(3) 



The last boundary condition is given by the temperature 
T(R^). It is obtained by a simultaneous solution of the 
temperature and density structure for the photosphere 
and the wind as described in the following. 

In calculating the atmospheric structure for the MS 
phases we closely follow de Koter (1993) and de Koter 
et al. (1995). Basically the atmosphere is characterised by 
two parts: the subsonic regime with an extended photo- 
sphere, and the wind, where the flow is accelerated to 
the terminal wind velocity v^. In between, both parts 
are smoothly connected. More precisely, in the subsonic 
regime the density structure p(r) is computed by solving 
the momentum equation for a stationary flow taking into 
account gas and radiation pressure. We write the momen- 
tum equation as 



dv 



1 dP g GM [1 - T(r) 



dr p dr 
where 
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(4) 



(5) 



is the ratio between the radiative acceleration gn and the 
gravitational acceleration g+ = GM/R 2 . The other vari- 
ables have their usual meaning. Note that T(r) is not con- 
stant due to its dependence on the flux weighted opacity 
Kp, which is a priori a function of T, p, v and dv / dr if the 
effects of continuum and line acceleration are taken into 
account. 

With regard to the flux mean opacity, we have chosen 
the following approach: For optical depths tr oss < t c = 2, 
we adopt the electron scattering opacity, i.e. up = <r e . 
Since <r e depends on the chemical composition and the 
ionization of the considered elements, we adopt the value 
calculated from the equation of state of the interior model 
at the boundary r*. At large optical depths (tr oss > t c ) we 
adopt for Kp(r) the Rosseland opacities as for the interior. 
Note that these opacity calculations assume LTE and a 
static medium. By checking the departure coefficients in 



depths LTE is attained to sufficient precision. Since opac- 
ities in expanding media are larger than in the static case 
(e.g. Karp et al. 1977), and the use of Rosseland opacities 
instead of flux weighted mean opacities can also lead to an 
underestimate of the opacity, the adopted opacities pro- 
vide a lower limit. An improved treatment would clearly 
require a much more complicated calculation of the radia- 
tive forces coupled together with a hydrodynamic solution. 

Following de Koter et al. (1995) we rewrite Eq. (4) us- 
ing the continuity equation and write the equation of state 
for the perfect gas as P g = a 2 p, where a = (kT / pinn) 1 ! 2 
is the local isothermal sound speed. This yields 



dv 
or 



1 



2a 2 



GM [1 - T(r 



k dT 



prriH or J 

The mean molecular weight p is obtained from the equa- 
tion of state as described for the electron scattering opac- 
ity above. Integrating Eq. (6) from subsonic velocities 
outward one progressively obtains an increasing velocity 
gradient, as v < a and the quantity in large brackets in 
Eq. (4) is negative for v < a. This property allows us to 
smoothly connect the subsonic part with a wind structure 
given by 



v(r) 



r 



(7) 



where the radius r is adjusted such that both the velocity 
field and its first derivative are continuous. The choice of 
the wind parameters and [3 is described in Sect. 2.4. 

A consistent solution of the momentum equation and 
temperature structure yields the density and temperature 
structure. At the same time the velocity is also ad- 
justed to fulfil the boundary condition with the stellar in- 
terior (cf. Eq. (1)). The temperature structure is given by 
radiative equilibrium in an extended grey atmosphere. It 
is determined from the generalised Eddington approxima- 
tion following Lucy (1971) and Wessolowski et al. (1988). 
The procedure is basically the same as described in SS94a, 
where more details can be found (cf. also de Koter 1993). 
In the outermost regions of the wind the temperature is 
not allowed to drop below a minimum value T m ; n , which, 
following the results of Drew (1989), is chosen as follows: 
For T e ff > 20 kK we set T mm = 0.4 T e ff, and T mm =0.5 T e ff 
for T eff < 20 kK. 

Because of the large computational effort required, we 
have neglected the influence of line blanketing on the tem- 
perature structure in calculating CoStar models. However, 
the effect of line blanketing is taken into account in the 
statistical equilibrium and radiative transfer calculations 
(see Sect. 2.5). 

2.3. Consistent interior and atmosphere solution 

For given stellar and wind parameters (L^, R+, M and M , 
i>oo, /?) we obtain with the procedure of Sect. 2.2 a consis- 
tent atmospheric structure covering the photosphere and 



ture T(it^) determining the temperature boundary condi- 
tion for the stellar interior. In another iterative process, 
embracing both the stellar interior and the atmosphere, we 
finally obtain a consistent solution for the entire star satis- 
fying all boundary conditions. The procedure is schemat- 
ically shown in Fig. 1. 

2.4. Mass loss rates and wind structure 

The adopted mass loss rate M evo i and the additional pa- 
rameters required to describe the wind structure are as 
follows: 1 . 

— Mass loss rates are adopted as in Meynet et al. (1994). 
This means that for population I stars throughout the 
HR diagram we use the mass loss rates given by de 
Jager et al. (1988), enhanced by a factor of two. Justifi- 
cations for this choice are given by Meynet et al. (1994) 
and Maeder & Meynet (1994). 

— The terminal velocities are from wind models of 
Leitherer et al. (1992). Comparisons of our adopted 
terminal velocities with observations are discussed in 
Sect. 3.2.1. 

— For the rate of acceleration of the supersonic flow (see 
Eq. (7)), we take [3 = 0.8 following theoretical pre- 
dictions of Friend & Abbott (1986) and Pauldrach 
et al. (1986). These predictions are in good agree- 
ment with observations of O stars by Groenewegen & 
Lamers (1991). 

2.5. non-LTE radiation transfer including line blanketing 

The non-LTE radiation transfer calculations, which yield 
the detailed spectral evolution use the atmospheric struc- 
ture from the CoStar model described above. For the de- 
tailed transfer calculations we used the ISA-WIND non- 
LTE code of de Koter et al. (1993, 1995). In this code, 
the line transfer problem is treated using the Sobolev ap- 
proximation, including the effects of the diffuse radiation 
field, and the continuous opacity inside the line resonance 
zone. As a new feature of the ISA-WIND code we also in- 
clude line blanketing, following the opacity sampling tech- 
nique introduced by Schmutz (1991). The method involves 
a Monte Carlo radiation transfer calculation including the 
most important spectral lines of all elements up to zinc. 
The ionization and excitation of the metals is treated as 
in SS94ab, where the reader is referred to for a detailed 
description of the entire procedure. 

The input physics for the atmospheric structure calcu- 
lations consists of atomic data for the elements explicitly 
included in the non-LTE model. In the present work hy- 
drogen and helium are treated using the same data as de 




Note that the mass loss rate M evo i adopted for the evolu- 
tionary calculations should not be confused with the theoreti- 
cal mass loss rate estimate M ca i c predicted by our atmosphere 
calculations in Sect. 3.2.2 



Fig. 2. a HR-diagram covering the MS phases for initial 
masses of 40, 60, and 85 Mq. The WR stage during the 
H-burning phase of the 85 Mq model is excluded. Crosses, 
triangles and stars denote the selected models on the 85, 60, 
and 40 Mq tracks respectively, for which detailed spectra have 
been calculated (see Paper II). The square denotes the position 
of an additional model presented in Paper II. b log (/-log T e jf 
diagram corresponding to the upper panel 



Koter et al. (1995): In the atomic model for H and He II 
we account for the first ten levels with principal quantum 
number n = 1 to 10. The atomic model for He I, consising 
of 17 levels, is described by Wessolowski et al. (1988). In 
total we account for 226 line transitions. The inclusion of 
C, N, O, Si, and other elements is in progress. 

The H, He, C, N, and O composition of the atmosphere 
is that corresponding to the outermost layer of the inte- 
rior model. For the metals included in the line blanketed 
atmosphere, the abundances of Anders & Grevesse (1989) 
have been adopted. 



3. Interior evolution and feedback to the ISM 

We have calculated three CoStar tracks for solar metallic- 
ity and initial masses of M; = 40, 60, and 85 Mq. Only 
the results for the MS phase are discussed in this work. 
We also exclude those parts of the tracks where the WR 
phase (defined by a hydrogen surface abundance X < 0.4 
in mass fraction and logT e ff > 4.; cf. Schaller et al. 1992) 
is already entered during the H-burning phase. In our 
models this occurs for the 85 Mq model. Therefore only 
part of the H-burning phase of this model is covered in 



Table 1. H-burnmg lifetimes and total momentum and energy deposition during MS evolution at Z=0.02 
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Fig. 3. a Mass loss rate as a function of age during the MS 
evolution. The line styles are as follows: 85 Mq short-dashed, 
60 Mq solid, and 40 Mq long-dashed, b Terminal velocity as 
a function of age. The lines are coded as in panel a 



the present publication. Later phases, including the Wolf- 
Rayet stages, are discussed in Schaerer (1995a, b). 

The evolutionary tracks of the models are presented 
in Fig. 2. Both the HR-diagram, and the corresponding 
gravity-effective temperature diagram are shown. The H- 
burning lifetimes are given in Table 1. As stated above, the 
Mi = 85 Mq model already enters the H rich WR phase 
during H-burning. However, for this mass the "normal" 
O star phase covers 90 % of the total MS phase (cf. Table 
1). A comparison of the lifetimes with the results from 
Meynet et al. (1994) (who use the usual plane parallel 
grey atmosphere as the outer boundary condition) shows 
that our 40 and 60 Mq models differ by less than 0.2 %. 
The largest difference, ~ 1.5 %, is obtained for the 85 Mq 
model. This shows that the overall MS evolution for these 
stars is not affected by the treatment of the atmosphere. In 



Section 4 we will examine in more details the atmospheric 
structures to understand this result. 

The evolution of the wind properties during the MS 
evolution are presented in Fig. 3, where we show the mass 
loss rate, M evo i, and the terminal velocity, i^, as a func- 
tion of age. Due to an increasing luminosity and evolu- 
tion towards lower temperatures M evo i increases during 
the evolution. The terminal velocity remains nearly con- 
stant during the first million years, whereas it decreases af- 
terwards, mainly as a result of an increasing radius, hence 
a decreasing escape velocity. 
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Fig. 5. Adopted terminal velocities as a function of es- 
cape velocity (The thick lines denote our models: line styles 
as in Fig. 3). Thin lines represent different fits. The thin 
dashed-dotted line shows the fit through our model data (see 
text). The dotted line (marked GLP) is the relation obtained 
by Groenewegen et al. (1989) from UV fits. See text for a dis- 
cussion 




Fig. 4. a Momentum flux of the 85 (short-dashed), 60 (solid), and 40 Mq (long-dashed) model as a function of its age. b 
Mechanical wind power emitted as a function of the age. The lines are coded as in the panel a 



3.1. Deposition of momentum and energy 

To study the effect of mass loss in massive stars on the in- 
terstellar medium it is of interest to quantify the wind mo- 
mentum flux Mevoi^oo , and the rate of release of mechani- 
cal energy 1/2 M^oxv^ into their surroundings. These val- 
ues may be used to investigate the properties of individ- 
ual nebula?, but may also be used to derive predictions 
for integrated young stellar populations (see Leitherer et 
al. 1992, Williams k Perry 1994) 

The momentum flux and mechanical wind power dur- 
ing the MS phase is shown in Fig. 4, while the total mo- 
mentum and energy deposition integrated over the MS 
lifetime is given in Table 1. Figure 4 illustrates the pro- 
gressive increase of both wind momentum and energy dur- 
ing the MS evolution caused by the strong enhancement 
of mass loss towards the end of the main sequence. 

3.2. Discussion of wind properties 

In this section, we discuss the most important adopted 
wind parameters, i.e. terminal velocity and mass loss rate. 
We first compare our adopted with observations. We 
then present theoretical mass loss rates derived from sim- 
ple energetic considerations, and compare these with the 
adopted values. 

3.2.1. Terminal velocities 

To make a comparison between the adopted and observed 
values of , we plot both against the escape velocity 
i> esc - The relation between i> esc and is discussed by e.g. 



Castor et al. (1975) and Abbott (1978). We derive the es- 
cape velocity from the effective gravity g e ff = g*{l — T). 
The correction for the radiation pressure due to electron 
scattering is given by T = 7.66 10 _5 <7 e (L/Lq) (M q /M) , 
where the value of the electron scattering opacity corre- 
sponds to the value at the boundary R+. (cf. Sect. 2.2). 

In Fig. 5 the adopted terminal velocity is shown as 
a function of the resulting i> esc . A least-square fit to 
the values from the three tracks shown in Fig. 5 yields 
Vco = (633.97 ± 4.42) + (1.921 ± 0.005) v esc , where the 
velocities are in km s _1 . As a comparison we have also 
plotted in Fig. 5 the relation derived by Groenewegen et 
al. (1989, Eq. 6) from UV fits including the effect of tur- 
bulence. The agreement is reasonable, the maximum dif- 
ferences in the considered range being about 20 %. Note 
that the overestimation obtained by their Eq. 6 for stars 
close to the ZAMS (i.e. large v esc ) would be slightly re- 
duced by adopting their fit-formula (Eq. 9), which takes 
also the dependence of on T e ff into account. For ve- 
locities 

^oo ^ 2000 km s our adopted formula probably 
slightly overestimates the wind velocity (cf. Leitherer et 
al. 1992). Our results are also in agreement with the rela- 
tion from Prinja et al. (1990), provided a small (~ 8 %) 
downward correction of their adopted escape velocities is 
applied. 

In resume, we can conclude that for escape velocities 
t'esc *J 800 km s _1 the adopted terminal velocities com- 
pare well with the observations, while for lower values the 
adopted is probably overestimated by up to ~ 25 %. 
However, as can be seen from Fig. 3, this only concerns a 
short period of time close to the end of the main sequence 
evolution. 
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Fig. 6. a Logarithm of the ratio of the theoretical mass loss rate M ca i c to the adopted value M as a function of effective 
temperature. Crosses, triangles and stars denote the selected models on the 85, 60, and 40 Mq tracks respectively. For a 
discussion see Sect. 3.2.2. b Theoretical wind momentum efficiency r] as a function of luminosity (same symbols as in panel 
(a)). Dashed line: Fit to the model predictions (Eq. (8)). The filled square indicates the result obtained from a consistent 
hydrodynamic model for £ Puppis (SS94a, model A). Dotted lined: Wind efficiency predicted by Lamers & Leitherer (f993, 
Eq. 23). For our main sequence models the predicted wind efficiencies are between ~ 0.2 and 0.7, which is up to a factor of 2 
larger than the values calculated by Lamers & Leitherer using the line force parameters of Pauldrach et al. (f 990) 



3.2.2. A mass loss rate estimate using the photon energy 
balance 

Presently all evolutionary models, including the ones de- 
veloped in this work, rely on the use of empirical mass loss 
rates. While radiation driven wind theory is quite success- 
ful in explaining the overall properties of OB stars and 
possibly also LBV's (e.g. Kudritzki et al. 1991 and refer- 
ences therein; Pauldrach & Puis 1990) discrepancies still 
remain (e.g. Groenewegen et al. 1989, Schmutz & Schaerer 
1992, Lamers & Leitherer 1993 hereafter LL93, Puis et 
al. 1995). Although the calculation of consistent hydrody- 
namic wind models, such as the ones presented by SS94ab, 
is clearly beyond the scope of the present work, it is how- 
ever very interesting to make estimates of the mass loss 
rates which can be driven by radiation pressure. 

We estimate the mass loss rate, subsequently called 
M ca \ c , from the photon energy balance, which is ob- 
tained from our Monte-Carlo radiation transfer calcula- 
tions. Since our models cover the entire main sequence for 
stars from M; = 40 to 85 Mq , a comparison of estimated 
radiation driven mass loss rates with observations could 
yield useful insight to understand the present difficulties 
of the radiation driven wind theory. 

To determine M ca \ c , we follow the considerations of 
Abbott & Lucy (1985). From the total radiative energy 
deposition in the wind, Lt = L^R*) — L(co), we calcu- 
late M ca i c from Lt = 1/2 Mcax^v^ + v^ sc ], assuming that 



the entire radiative energy deposition in the wind is con- 
verted to mechanical energy only, by lifting matter out of 
the gravitational field and giving it its asymptotic kinetic 
energy. The flux transfer rate Lt is obtained from the MC 
simulation taking into account a large number of metal 
lines, and including line overlap and multiple scattering 
(see SS94a and Paper II). Note, however, that achieving 
consistency is beyond the scope of the present work. This 
could be accomplished by adjusting the mass loss rate un- 
til Lt, which depends on the wind density, equals the 
mechanical energy (cf. Abbott & Lucy 1985). Our results 
should only therefore be taken as estimates for mass loss 
driven by radiation pressure. Surprisingly, as shown below, 
this method shows a good agreement with detailed hydro- 
dynamic calculations of Schaerer & Schmutz (1994a). 

In Figure 6a, we have plotted the ratio M ca \ c / M evo \ of 
the theoretical to the adopted mass loss rate as a func- 
tion of effective temperature. For the models close to the 
ZAMS (T e flf ^ 43.6 kK) we see that the estimated theoret- 
ical mass loss rate is of the same order, or even larger, than 
the adopted values for M evo i. In this case, the energy ex- 
tracted from the radiation field can, in principle, account 
for the energy stored in the stellar wind. On the other 
hand, for the models with T e ff 43.6 kK, the adopted 
M evo i is larger than the theoretical derived values. Sim- 
ilar discrepancies as a function of effective temperature 
are predicted by LL93, who use analytic solutions of Ku- 



Pauldrach et al. (1990). This indicates a shortcoming in 
the current state of the radiation driven wind theory for 
evolved stars, which may be resolved by consistent hydro- 
dynamic calculations. 

Of particular interest is the wind efficiency rj = 
Mvoo /(L^/c) predicted by radiation driven wind models. 
For OB stars one typically finds rj ~ 0.1-0.6, while for 
WR stars values considerably larger than unity can be 
obtained (e.g. LL93). 

In Figure 6b, we have plotted the theoretical wind ef- 
ficiency ?y ca i c as a function of the luminosity. The model 
data are fitted by the relation 

log^aic = (0.669 ± 0.067) log(i/^0) - 4.276 ± 0.386, (8) 

with a rms of 0.052 dex. The relation is represented by the 
dashed line. Interestingly the corresponding value from a 
self-consistent hydrodynamic calculation of SS94a for £ 
Puppis is quite well matched by the above relation. This 
indicates that, despite the lack of consistency pointed out 
above, the estimated mass loss rate M ca \ c determined from 
the flux transfer rate Lt yields a reasonable value for the 
mass loss rate as derived from full hydrodynamic model- 
ing. 

We have compared our predictions for ?y ca i c with the 
those of LL93, which are based on parametrised line forces 
calculated by Pauldrach et al. (1990). The results of LL93 
(their Eq. 23) are plotted as the dotted line in Fig. 6b, 
the uncertainty being 0.1 dex. Figure 6b reveals that our 
models predict a steeper increase of the wind momentum 
with luminosity. For the most luminous model we obtain 
??caic ~ 0.7, which is a factor of two larger than the value 
of LL93. Since the methods used in both approaches are 
completely independent, it is difficult to trace the differ- 
ences back to one single reason. However, we expect to 
find the largest differences for cases where the effects of 
multiple scattering may become important since this effect 
is correctly treated in our models, but has been neglected 
by Pauldrach et al. (1990). On the other hand the larger 
wind efficiencies predicted by our calculations could also 
indicate systematic differences in the ionization structures. 
With respect to the ionization problem, we note that our 
models are supported by comparisons of predicted Fe fea- 
tures in UV spectra, which show a good agreement with 
observations for the evolved models (see Paper II), where 
iron is the dominant source of the radiation force. 

New calculations of radiation driven wind mod- 
els based on the improvements reported by Pauldrach 
et al. (1994) have recently been presented by Puis et 
al. (1995). They show that the most appropriate quantity 
when comparing theoretical predictions with observations 

1/2 

is the "wind momentum rate" Mv^R^ , since this quan- 
tity is expected to show a very weak dependence on the 
adopted stellar parameters. Indeed a strong correlation of 
log ( MvqqR]/ 2 J with log L is expected if the line force 
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Fig. 7. Theoretical wind momentum rate as a function of lumi- 
nosity. Stars, triangles and crosses denote the predictions from 
our models on the 40, 60, and 85 Mq tracks respectively. Small 
squares show the predictions for individual galactic objects of 
luminosity class I (filled squares) and II-V (open squares) from 
Puis et al. (1995). The lines show the mean relations for ob- 
served Luminosity class I (solid) and II-V objects (dashed) 
derived from the data of Puis et al. (their Table 8, galactic 
objects only) 



parameter k, which represents the flux- weighted number 
of driving lines, is constant. The slope is then only deter- 
mined by the line force parameters a and 8, and is found 
to be l/(a — 8). Observational evidence for such a corre- 
lation has been presented by Kudritzki et al. (1995). 

To compare our calculations with the results of Puis et 
al. (1995) we have plotted the predicted wind momentum 
rate from our models as a function of luminosity in Fig. 7. 
Also shown are the predictions for the individual galac- 
tic objects from Puis et al. and the mean observed values 
for luminosity class I objects (solid) and the objects of 
class II-V (dashed). Figure 7 shows that our models re- 
produce well the observed behaviour, i.e. both the slope 
of the relation and the trend of increased wind momen- 
tum rate for supergiants. On the average our predictions 
seem to be in rough agreement with the results of Puis 
et al. (1995). However, as pointed out by Puis et al. and 
apparent in Fig. 7, their models show a discrepancy for 
the supergiants as their theoretically predicted momen- 
tum is about 0.25 dex too small. This discrepancy is not 
present in our models, hence our supergiant models are in 
better agreement with observations. Whether this differ- 
ence with the recent Munich models is due to the effect 
of multiple scattering, which in contrast to Puis et al. is 
included in our calculations, or has other causes (different 



careful investigation. As mentioned above iron, which is 
the dominant source of the radiation force, is, however, 
well described in our supergiant models (see Paper II). 
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Fig. 8. HR-diagram for the initially 60 Mq model, calculated 
with CoStar (solid line) and the standard Geneva evolution 
code (dashed). In both cases the boundary with the stellar 
interior model is set at rR O ss=20. The standard track (plane 
parallel atmosphere, tross = 2/3) is shown with the dotted 
line. Also plotted is the zero-age main sequence (ZAMS). The 
difference between the tracks is negligible in view of the errors 
related to observational determinations of temperatures and 
luminosities 



4. Comparison of plane parallel and extended at- 
mospheres for evolutionary models 

In this section we discuss the influence of atmospheric 
boundary conditions on the evolutionary tracks. As will 
be discussed below, we need to adopt a slightly different 
treatment for the most massive stars (M; ^ 85M ), which 
lie close to the Eddington limit. We therefore discuss these 
cases separately. 

4-1. Evolution of the Mi = 40 & 60 Mq models 

In Fig. 8 we show a detailed comparison between the 60 
Mq track calculated with the conventional atmospheric 
treatment and the stellar models including spherically ex- 
panding atmosphere ( Co Star models). To compare exactly 
the same definitions of effective temperatures, and hence 
radii, we have set the boundary at tr oss =20 for both 
cases. Different tracks should thus reflect the influence 
of different boundary conditions between plane parallel 



lated with the standard treatment, i.e. plane parallel at- 
mosphere with the inner boundary at tr oss = 2/3 (dotted 
line). 

As we can see from Fig. 8 the tracks hardly differ. 
Similarly, the internal structures are basically identical, 
which explains why the H-burning lifetimes do not differ. 
This is due to the fact that the wind is mostly optically 
thin, which implies that the thermal conditions are es- 
sentially determined in the quasi-hydrostatic photosphere 
and hence the boundary conditions including the wind are 
nearly identical to those of plane parallel models. 

With respect to the standard track, both the plane 
parallel and spherically expanding models with the inner 
boundary set at large optical depth (tr oss =20) evolve at 
slightly higher temperature during most of the MS. The 
difference amounts to only AT <> 400 K, if the same ages 
are compared. On the other hand the spherical extension, 
measured by the ratio of radius R2/3 at tr oss = 2/3 to 
reaches only up to R2/3/R* 1.015 for the most evolved 
models. This confirms the expectation that despite the 
spherical extension of the atmosphere the radius definition 
causes no ambiguity for the MS phase of stars with M; = 
40&60M Q . 

4.2. Evolution of the Mi = 85 Mq model 

For more massive stars the situation can become slightly 
more complicated because they evolve quite close to the 
Eddington limit. This may introduce an additional diffi- 
culty in the computation of the atmospheric structure. 

To illustrate this point we have calculated the MS 
evolution of a 85 Mq model with both types of atmo- 
spheres. We first assumed a constant mass loss rate of 
M evo i = 1.14 10 -5 MQyr -1 to eliminate indirect effects of 
mass loss on the evolutionary tracks (see Sect. 4.2.2). The 
resulting tracks are shown in Fig. 9a. It is important to 
note that, contrary to the usual treatment (cf. Sect. 2.2), 
the atmosphere of the CoStar model has been calculated 
with a depth-independent effective gravity corrected for 
electron scattering only. As we will argue below, this 
should yield a more realistic structure than the hydrostatic 
plane parallel atmosphere using Rosseland opacities. 

4.2.1. Constant mass loss tracks for Mi = 85 Mq 

Figure 9b shows the evolutionary track in the gravity- 
Teg- diagram. To illustrate its proximity with respect to 
the Eddington limit we have plotted the limit derived 
by Lamers & Fitzpatrick (1988) from Kurucz models (la- 
beled <7Edd)- Also shown is the lowest log 3 value (labeled 
log <7min ) for which they obtained converged hydrostatic 
atmosphere structures, and which were used to determine 
the <7Edd limit by extrapolation. 

The difference between the tracks in the HR-diagram 
(Fig. 9a) can be understood by looking at the atmospheric 
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Fig. 9. a HR-diagram of a 85 Mq model calculated with constant M e voi = l-14 10 -5 MQyr -1 . The solid line shows the CoStar 
model adopting a depth-independent effective gravity. The dashed track is for the conventional plane parallel atmosphere. The 
CoStar track evolves at slightly higher T e jf due to a larger pressure at the boundary r*. The star marks a model discussed in 
detail (cf. Fig. fO). b log g vs. log T e jf diagram for the 85 Mq track calculated with constant M evo i (solid line). The star marks 
the position of the same model as in panel (a). As a comparison the track of the 40 Mq model from Fig. 2 is also plotted (long 
dashed). The 60 Mgtrack is not shown here, since it virtually coincides with the dotted line. The short dashed line, marked 
</Edd, traces the Eddington limit determined by Lamers & Fitzpatrick (1988) from Kurucz models by extrapolation from the 
lowest gravity models shown by the dotted line (c/ m in). This illustrates the proximity of the 85 Mq track to the Eddington limit 



structures of models with stellar parameters correspond- 
ing to the model marked 306 in this figure. Figure 10 
shows the temperature and density stratification of the 
plane parallel (dashed line) and the spherically expanding 
atmosphere (solid line). Note that at log tr oss 2; —0.2 the 
hydrostatic model shows a small density inversion. The p— 
inversion occurs in the zone where He becomes completely 
ionised and is due to the increasing opacity (Fig. 9b). 

If we suppress the density inversion in the plane par- 
allel model by assuming a constant opacity (k = a e ), the 
density scale height remains constant (dotted line). How- 
ever, a higher density results at any given optical depth, 
and at the inner boundary in particular, because the opac- 
ity only accounts for electron scattering. 

In a hydrodynamic solution, on the other hand, one 
could make the conjecture that the opacity increase at 
rather low optical depths contributes to the acceleration of 
the outflow and washes out the density inversion. Since in 
this work we are not able to solve consistently for the full 
hydrodynamic equations of the entire atmospheric struc- 
ture, we use our usual procedure to describe the wind and 
the photosphere, but we do not allow for density inversions 
in CoStar models. This is simply obtained by adopting 
a constant, i.e. depth-independent effective gravity cor- 
rected for electron scattering only. The usual Rosseland 
opacities are however used for the temperature determi- 



nation. The resulting structure of the spherically expand- 
ing atmosphere is plotted in Fig. 10 showing (in the inner 
part) the same density scale height as the plane parallel 
model where the p-inversion has been suppressed. One 
can clearly see the large velocity gradient, which is just 
located outward of the zone of the density inversion of the 
hydrostatic model. 

Since the temperature is basically determined in the 
photosphere (cf. above) both models have the same tem- 
perature at the inner boundary. Due to the higher density 
in atmospheres without p-inversion the pressure is, how- 
ever, slightly larger 2 than in the hydrostatic atmosphere 
where the density inversion occurs. 

When the modified boundary conditions are used to 
determine the entire stellar structure, the star will read- 
just itself to the increased compression, which yields a 
slightly lower radius. This explains the differences be- 
tween the CoStar and the conventional track (both with 
M evo i=const) shown in Fig. 9. As we will see below, the 
effect of modified boundary conditions has indirect con- 
sequences on the evolution, since the position in the HR- 
diagram determines the mass loss rate. 



2 Approximately 30 % of the total pressure at is provided 
by the gas pressure. 





Fig. 10. Atmosphere structure for model 306 (cf. Fig. 9). Dashed line: plane parallel atmosphere; dotted line: also plane 
parallel, but adopting a depth-independent c/ e ff using k = a e ; solid line: spherically expanding atmosphere with constant c/ e ff in 
the photosphere, and M eV oi = l-14 10 -5 MQyr -1 . a Density structure as a function of Rosseland optical depth, b Temperature 
structure, c Rosseland opacity as a function of Rosseland optical depth, d Velocity gradient of the spherically expanding 
atmosphere as a function of optical depth. 



4.2.2. Variable mass loss track for M; = 85 Mq 

We now adopt the mass loss rate prescription given in 
Sect. 2.4. This allows us to illustrate indirect effects of 
the boundary conditions on evolutionary tracks. Figure 11 
shows a comparison of the CoStar track with the track cal- 
culated with the plane parallel atmosphere (from Meynet 
et al. 1994). Here we concentrate on the evolution from 
the ZAMS to the beginning of the WR phase, marked by 
a star in the HR-diagram. To illustrate the differences in 
later phases, we show the continuation of the tracks in the 
WNL phase. As explained above, the CoStar track has 
been calculated with a constant effective gravity in the 
photosphere. Therefore the pressure at the boundary 
is slightly larger than for the plane parallel atmosphere, 
which results in a reduced radius, i.e. a larger effective 
temperature. Since, in turn, the mass loss rate at a given 
luminosity decreases with increasing T e ff (see de Jager et 
al. 1988; also Lamers & Leitherer 1993), the CoStar model 
suffers a smaller mass loss (cf. Fig. lib), which explains 
its slight overluminosity and hence the shortened lifetime 
(although only ~ 1.5 %) with respect to the conventional 
model. 

Clearly the effects discussed here for the O star phase 
are quite small. Looking at the beginning of the subse- 
quent WR phase, part of which is also plotted in Fig. 11, 
the tracks seem to diverge. In fact, this only concerns a 
very short time before the end of H-burning and the be- 
ginning of the hydrogen free WR phase (WNE). 



We conclude that the differences in evolutionary track, 
interior evolution and lifetimes on the main sequence, be- 
tween the conventional and CoStar models are negligible. 
However, a small uncertainty in the tracks and the pre- 
dicted stellar parameters remains for the most luminous 
MS stars due to their proximity to the Eddington limit. 
A fully hydrodynamic treatment including the subphoto- 
spheric layers and the wind is required to improve the 
present treatment. 

Post-MS phases are discussed in Schaerer (1995a, b). 

5. Summary and conclusions 

In the present paper we have presented the first "com- 
bined stellar structure and atmosphere models" (CoStar) 
for massive stars, which consistently treat the stellar in- 
terior and a spherically expanding atmosphere including 
the wind. Our approach replaces the widely used boundary 
conditions given by plane parallel grey atmospheres. The 
CoStar models also predict the detailed emergent spec- 
trum along the evolutionary tracks taking non-LTE effects 
and line blanketing into account (see Paper II). 

As a first exploration of the behaviour of radiation 
driven wind models on the entire MS we derive esti- 
mates of theoretical mass loss rates from energy consid- 
erations (Sect. 3.2.2). This is of particular interest since 
our atmosphere calculations include the effects of mul- 
tiple scattering and line overlap, which are usually ne- 
glected. We have compared our results with the predic- 
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Fig. 11. a HR-diagram of the 85 Mq CoStar model (solid line) 
and the conventional model from Meynet et al. (1994; dashed 
line: uncorrected temperature T*). The H-burning phase pro- 
ceeds from the ZAMS up to the position marked with a star 
where the tracks enter the WR phase (WNL). Since the CoStar 
model evolves at slightly higher temperatures its mass loss is 
smaller (cf. panel (b)) and hence it evolves at higher lumi- 
nosities, (b): Evolution of the central hydrogen abundance X c 
(mass fraction) and the stellar mass, expressed by the ratio 
M I 'Mi of the present mass to the initial mass. Shown is only 
the evolution from the ZAMS up to the beginning of the WR 
phase, which is marked by a star in the panel (a). The mass 
difference at this stage is PS 2.8 Mq. Although the tracks di- 
verge during the WR phase, the resulting differences for the 
H-burning lifetimes are still small 



tions from recent wind models of the Munich group (Paul- 
drach et al. 1990, 1994, Puis et al. 1995). While we find 
an overall agreement with their results, our models in ad- 
dition also reproduce the strong wind momentum rates 
observed in supergiants. Possible reasons for this finding 
have been discussed. Consistent hydrodynamic calcula- 
tions (cf. Schaerer & Schmutz 1994a) will be necessary to 
verify whether for OB stars the quantitative discrepancies 
of the wind theory (see e.g. Puis et al. 1995) can indeed 
be explained by the effects included in our models. 

One of the important aims of this study was to inves- 
tigate the possible influence of the spherically extended 
atmosphere on the main sequence evolution of massive 
stars, and in particular, on the predicted positions in the 
HR-diagram. From the well known convergence proper- 
ties of radiative envelopes (e.g. Schwarzschild 1958) it is 
expected that small changes of the external boundary con- 
ditions have little influence on the interior evolution dur- 
ing phases where only a central nuclear burning source 
is present. This is confirmed in general by our model- 
ing of the MS evolution, although some uncertainties, re- 
lated to the modeling of either density inversions in plane 
parallel hydrostatic atmospheres or the acceleration zone 



(Sect. 4). 

It is important to realize that the above stated rela- 
tively unimportant effect of outer boundary conditions on 
the stellar structure may not generally be true in phases 
where a burning shell is present. In this case the bound- 
ary values may indeed influence the conditions in the shell, 
leading to a readjustment of its nuclear energy production, 
which thereby also affects the growth of the underlying 
He-core. Such a behaviour was e.g. found for stars with 
M ~ 15-20 Mq undergoing blue loops (Langer 1991). 
Subsequent to the MS evolution presented in this paper 
the considered models will evolve through a brief "LBV 
like" phase with strong mass loss before entering the Wolf- 
Rayet phase (see Meynet et al. 1994), where the burning 
shell will be extinguished in most cases. While the effects 
of the expanding atmosphere on the interior evolution has 
been studied for the WR phases by Schaerer (1995a, b), its 
influence during the short transient phase is expected to 
be negligible, although it still needs to be explored. 

With respect to the accuracy of positions in the HR- 
diagram the situation can be summarised as follows: Since 
for main sequence OB stars the continuum is essentially 
formed in a quasi-hydrostatic photosphere, not only the 
subsonic lower boundary conditions of both spherical and 
plane parallel atmospheres are basically identical, but also 
the spherical extension is negligible. The predicted stellar 
parameters (radius, T e ff) are therefore essentially identical 
(Sect. 4). 
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